
# Replication materials ---------------------------------------------------
# Paper:   Estimating Uncertainty in Party Policy Positions 
#          Using the Confrontational Approach
# Authors: Tom Louwerse, Huib Pellikaan
# Contact: Tom.Louwerse@tcd.ie
# Dependencies: JAGS, script 'confrontational.R', 'confrontational_extra_
#               functions.R', packages rjags, foreach, plotrix

# Instructions for use ----------------------------------------------------
# 1. Make sure you have JAGS installed on your system
#    Download latest version from: 
#    http://sourceforge.net/projects/mcmc-jags/files/JAGS/
# 2. Install packages rjags, foreach and plotrix:
#    install.packages(c("rjags","foreach", "plotrix"))
# 3. Set working directory. All data and the 'confrontational.R' script
#    should be available in that directory.


# Working directory -------------------------------------------------------
# setwd("")

# Load packages -----------------------------------------------------------
library(foreach);library(plotrix);library(rjags)

# Source confrontational function -----------------------------------------
source("confrontational.R")
source("confrontational_extra_functions.R")

# Load data ---------------------------------------------------------------
dta_nl <- local(get(load("Data Confrontational NL 2010.RData")))
load("NL2010_comparison_data.RData")
load("Data EU.RData")


# Descriptives ------------------------------------------------------------
descriptive_data <- foreach(i=1:length(dta_nl), .combine=rbind) %do% {
  data.frame(dimension=names(dta_nl)[i],n.parties=nrow(dta_nl[[i]]),n.items=ncol(dta_nl[[i]]), contra=table(dta_nl[[i]])[1], neutral=table(dta_nl[[i]])[2], pro=table(dta_nl[[i]])[3])
}

write.csv(descriptive_data, file="descriptives.csv")

# IRT model (NL) ----------------------------------------------------------
model_nl = foreach(i=1:length(dta_nl)) %do% confrontational(dta_nl[[i]], seed=5)

# Create plot (Figure 1)
windows(9,11)
par(mfrow=c(4,3), mar=c(2, 2, 3.5, 2))
for(i in 1:length(model_nl)) {
  plot(model_nl[[i]], main=names(dta_nl[i]))
} 


# Diagnostics -------------------------------------------------------------
foreach(i=1:length(model_nl)) %do% heidel.diag(model_nl[[i]]$actors)
foreach(i=1:length(model_nl)) %do% raftery.diag(model_nl[[i]]$actors)
foreach(i=1:length(model_nl)) %do% geweke.diag(model_nl[[i]]$actors)



# Comparison with CMP CIs -------------------------------------------------
# Select only dimensions on which we have priors from CHES
dta.set.select <- match(names(dta.ConCMP.per.issue), names(dta_nl))

# Get correct order of parties in original dataset
dta.set.order <- match(dta.ConCMP.per.issue[[1]]$Party, rownames(dta_nl[[1]]))

# Parse CMP result into a dataframe, so it can be plotted using plot.confrontational
dta.CMP.result <- foreach(i=1:length(dta.ConCMP.per.issue)) %do%
  list(actors_summary=data.frame(mean=dta.ConCMP.per.issue[[i]]$CMP, 
             low=dta.ConCMP.per.issue[[i]]$CMP-(1.96 * dta.ConCMP.per.issue[[i]]$CMP_se),
             hi=dta.ConCMP.per.issue[[i]]$CMP+(1.96 * dta.ConCMP.per.issue[[i]]$CMP_se),
             row.names=dta.ConCMP.per.issue[[i]]$Party))

# Create plot (Figure 2)
windows(7, 10)
par(mfrow=c(3,2), mar=c(2, 2, 3, 2))
CMP.dimensions <- c("Multiculturalism", "Stateconomy", "Traditional morality", "Liberalism-Conservatism", "Nationalism", 
                    "Liberalism-Conservatism", "EU", "Environment", "Militarism", "State Services")
for(i in c(1,2,3,7,10)) {
  biplot(model_nl[[dta.set.select[i]]], dta.CMP.result[[i]], labelx="Confrontational Approach", labely="Manifesto Project", main=paste(names(dta_nl)[dta.set.select[i]], CMP.dimensions[i], sep=" / "))
}



# IRT model with priors ---------------------------------------------------

# Select only dimensions on which we have priors from CHES
dta.set.select <- match(names(dta.ConCHES.per.issue), names(dta_nl))

# Get correct order of parties in original dataset
dta.set.order <- match(dta.ConCHES.per.issue[[1]]$Party, rownames(dta_nl[[1]]))

# Run those models
model_nl_priors = foreach(i=1:length(dta.set.select)) %do% 
  confrontational(dta_nl[[dta.set.select[i]]][dta.set.order,],
                  priors=list(mean=dta.ConCHES.per.issue[[i]]$CHES,
                              se=dta.ConCHES.per.issue[[i]]$CHES_se),
                  seed=5)


# Compare CI width (models with and without priors) (reported in-text)
calcCIwidth <- function(x,y) {
  p <- merge(x,y,by=0)
  CI.width.normal <- p$'97.5%.x' - p$'2.5%.x'
  CI.width.prior <-  p$'97.5%.y' - p$'2.5%.y'
  return(data.frame(Party=p$Row.names,
                    CI.width.normal=CI.width.normal,
                    CI.width.prior=CI.width.prior,
                    Difference=CI.width.normal-CI.width.prior))
}

ci.width <- foreach(i=1:length(model_nl_priors)) %do% 
  calcCIwidth(model_nl[[dta.set.select[i]]]$actors_summary, model_nl_priors[[i]]$actors_summary)
names(ci.width) <- names(dta_nl)[dta.set.select]

ci.width.diff.av <- foreach(i=1:length(ci.width), .combine=rbind) %do% 
  data.frame(Issue=names(ci.width)[i],
             Mean_Difference=mean(ci.width[[i]]$Difference),
             Perc_Difference=mean(ci.width[[i]]$Difference)/mean(ci.width[[i]]$CI.width.normal) * 100)

print( ci.width.diff.av)

# Correlations between models with and without priors (reported in-text)

corModels <- function(x,y,method="pearson") {
  p <- merge(x,y,by=0)
  return(cor(p$x,p$y, method=method))
}

ci.cor.r <- foreach(i=1:length(model_nl_priors), .combine=c) %do% corModels(model_nl[[dta.set.select[i]]]$actors_summary[,1], model_nl_priors[[i]]$actors_summary[,1])
names(ci.cor.r) <- names(dta_nl)[dta.set.select]

ci.cor.tau <- foreach(i=1:length(model_nl_priors), .combine=c) %do% corModels(model_nl[[dta.set.select[i]]]$actors_summary[,1], model_nl_priors[[i]]$actors_summary[,1],method="kendall")
names(ci.cor.tau) <- names(dta_nl)[dta.set.select]

print(ci.cor.r)
print(ci.cor.tau)



# Different specifications for null-scores --------------------------------

# Create dataset with zero values as missing (NA)
cul <- dta_nl[[1]]
cul[which(cul==2, arr.ind=TRUE)] <- NA

out_org <- confrontational(dta_nl[[1]],seed=5)
out_sum <- (rowSums(dta_nl[[1]])-20)
out_mar <- confrontational(cul,seed=5)
out_leg_cut <- confrontational(dta_nl[[1]], model.file="jags_model_actor_cut.bug",seed=5)

# Use function compareplot.confrontational
compareplot.confrontational(out_org, out_mar, out_leg_cut)


# IRT Model (EU) ----------------------------------------------------------

# Run model
model_eu <- confrontational(eu_dta, seed=5)

# Create plot (Figure 4)
windows(7,7)
par(mfcol=c(1,1),  mar=c(2, 2, 3, 2))
plot(model_eu)

# Scores per party group (not reported)
partygroups = list()
partygroupscores = matrix(NA,ncol=3,nrow=length(levels(eu_metadata$PARTYGROUP)))
partygroupnames = levels(eu_metadata$PARTYGROUP)
rownames(partygroupscores) = partygroupnames
colnames(partygroupscores) = c("Mean", "2.5%", "97.5%")
for(i in 1:length(levels(eu_metadata$PARTYGROUP))) {
  partygroups[[i]] = which(eu_metadata$PARTYGROUP==levels(eu_metadata$PARTYGROUP)[i])
  sc = apply(model_eu$actors, 1, function(x) mean(x[partygroups[[i]]]))
  partygroupscores[i,] = c(mean(sc), quantile(sc, 0.025), quantile(sc, 0.975))
}
print(partygroupscores)

# Scores per country (means and standard deviations)
france = which(eu_metadata$COUNTRY=="FR")
ireland = which(eu_metadata$COUNTRY=="IR")
netherlands = which(eu_metadata$COUNTRY=="NL")
countryscores = apply(model_eu$actors, 1, function(x) 
  c(weighted.mean(x[france],eu_metadata$SEATS[france]),
    weighted.mean(x[netherlands],eu_metadata$SEATS[netherlands]),
    weighted.mean(x[ireland],eu_metadata$SEATS[ireland])))
countryscores_sd = apply(model_eu$actors, 1, function(x) 
  c(sd(x[france]),sd(x[netherlands]),sd(x[ireland])))

# Calculate summary statistics, save in dataframe ready for plotting
country_means <- list()
country_means$actors_summary = t(apply(countryscores, 1, function(x) c(mean(x), quantile(x, 0.025), quantile(x, 0.975))))
rownames(country_means$actors_summary) = c("France", "Netherlands", "Ireland")

country_sds <- list()
country_sds$actors_summary = t(apply(countryscores_sd, 1, function(x) c(mean(x), quantile(x, 0.025), quantile(x, 0.975))))
rownames(country_sds$actors_summary) = c("France", "Netherlands", "Ireland")

windows(7,7)
# Plot (Figure 4)
par(mfcol=c(2,1),mar=c(2.1,4.1,2.1,2.1))
plot.confrontational(country_means,labelcex=1)
title(ylab="Mean")
plot.confrontational(country_sds,labelcex=1)
title(ylab="Standard deviation")


# Calculate posterior probabilities that country score row is higher than country score column (reported in-text)
countryscores_matrix = matrix(NA, nrow(countryscores),nrow(countryscores))
countryscores_sd_matrix = matrix(NA, nrow(countryscores),nrow(countryscores))

for(i in 1:nrow(countryscores)) {
  for(j in 1:nrow(countryscores)) {
    countryscores_matrix[i,j] <- sum(countryscores[i,] < countryscores[j,])/length(countryscores[i,])
    countryscores_sd_matrix[i,j] <- sum(countryscores_sd[i,] < countryscores_sd[j,])/length(countryscores_sd[i,])
  }
}
rownames(countryscores_matrix) = colnames(countryscores_matrix) = c("France", "Netherlands", "Ireland") 
rownames(countryscores_sd_matrix) = colnames(countryscores_sd_matrix) = c("France", "Netherlands", "Ireland") 

print(countryscores_matrix)
print(countryscores_sd_matrix)

